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In recent years deep X-ray and infrared surveys have provided an efficient way to find accreting 
supermassive black holes, otherwise known as active galactic nuclei (AGN), in the young uni- 
verse. Such surveys can, unlike optical surveys, find AGN obscured by high column densities of 
gas and dust. In those cases, deep optical data show only the host galaxy, which can then be 
studied in greater detail than in unobscured AGN. Some years ago the hard spectrum of the X- 
ray "background" suggested that most AGN were obscured. Now GOODS, MUSYC, COSMOS 
and other surveys have confirmed this picture and given important quantitative constraints on 
AGN demographics. Specifically, we show that most AGN are obscured at all redshifts and the 
amount of obscuration depends on both luminosity and redshift, at least out to redshift z ~ 2, 
the epoch of substantial black holes and galaxy growth. Larger-area deep infrared and hard 
X-ray surveys will be needed to reach higher redshifts and to probe fully the co-evolution of 
galaxies and black holes. 



1. Cosmic Growth of Black Holes and Galaxies 

Abundant evidence indicates that the growth of a supermassive black hole is closely 
tied to the formation and evolution of the surrounding galaxy. The energy released from 
accretion onto the black hole affects star formation in the galaxy, probably limiting 
growth at the high- and low-mass ends, and of course the distribution and angular mo- 
ment um of matter in th e galaxy governs the amount of m atter accumulated by the black 
hole ( Silk h Reesl Il998t iKind l2005t iRovilos et aLl l2007h . Emergent energy fr om accre- 
tion i s also a factor in understanding ionization and radiation backgrounds (jHasinger 
2000; lLawrencdl2001l ) . Understanding the growth history of these black holes is therefore 
critical to understanding the global evolution of structure in the Universe. 

Yet the demographics of supermassive black holes remain elusive. The largest samples 
of quasars and Active Galactic Nuclei (AGN), by which we mean supermassive black 
holes in a high accretion-rate phas^j], ha ve been found through optic al selection (e.g., the 
Sloan Digital Sky Survey quasar sample; ISchneider eFaZll2002l [2007h . but, at least locally, 
these are not representative of the larger AGN population. Instead we need surveys less 
biased against obscured AGN. 

There are three reasons to suspect that most AGN are obscured by large column densi- 
ties of gas and dust. First, a large body of evidence suggests that local AGN have geome- 
tries that are not spherically symmetric, and that different aspect angles present markedly 
different observed characteristics; this is referred to as AGN unification (IAntonuccilll993t 



Urrv fc Padovanilll995f ). Second, AGN are more common at high redshift (z ~ 2 — 3), 
where the average star formation rate is higher and thus it is even more likely that 
gas and dust surround the galaxy nucleus than at z ~ 0. Third, and most important, 
obscured AGN are required to explain the shape of the X-ray "background" radiation. 
The X-ray "background" is actually the superposition of individual AGN that were 

| Some make a distinction between AGN and quasars, with the latter being above an arbitrary 
luminosity, typically Mb — —23 mag. Here we use the term AGN to refer to an actively accreting 
supermassive black hole above or below this luminosity. 
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Figure 1. (Left) The X-ray background spectrum (shown in units of energy per logarithmic 
band) is very hard, peaking near 30 keV, in contr ast to the spectra of unobscured AGN which 
are roughly horizontal in these units. (Figure from Gilli 2004.) (Right) Large column densities of 
gas absorb low-energy photons via the photoelectric effect, such that emission from an obscured 
AGN peaks at an energy that increases with column density. Each line represents a factor 
of 10 increase in equivalent Nh (lin es are marked with log Nh in atoms/cm 2 ; assumes solar 
abundances and the cross sections 

of|M 

ornson fc McCammonlll983l ). 



not resolved in early X-ray experiments (hence the designation "background"). As shown 
in Figure la, its spectrum peaks at ~ 30 keV (i.e., this is where most of the energy is 
produced) , much harder than the typical spectrum (roughly flat in these units) of an un- 
obscured AGN. In obscured AGN, however, the softest X-ray photons have been absorbed 
via the photoelectric effect, and Compton-thick AGN (those with Nh ^ 10 24 cm~ 2 ) actu- 
ally peak at roughly 30 keV (Fig.QJ)). Thus X-ray observations have long indicated there 



Comastri et al 



is a l arge population of heavily obscured AGN (|Setti &; Woltierl 11989c 
1 19951: iGilli et al. 2001 ). While much of the X-ray background has been resolved at low 
energies, recent work on X-ray deep fields suggests the h ardest (most obscured) X-ray 



sources have yet to be detected (e.g.. IWorslev et aLl l2005) 



Relatively unbiased AGN samples require joint hard X-ray and infrared surveys. Hard 
X-rays (E > 2 keV) penetrate all but the thickest column densities and efficiently lo- 
cate black hole accretion. The absorbed optical through soft X-ray emission heats the 
surrounding dust and is re-radiated in the infrared, at wavelengths that depend on the 
dust temperature. Meanwhile, high-resolution optical surveys (e.g., with HST) allow sep- 
aration of nuclear (AGN) and host galaxy light. Today, NASA's three operating Great 
Observatories — the Chandra X-ray Observatory, the Spitzer Space Telescope, and the 
Hubble Space Telescope — enable matched X-ray, infrared, and optical surveys at un- 
precedented depth and resolution. 



2. Deep Multiwavelength Surveys 

The announcement in spring 2000 of the first Spitzer Legacy opportunity led to the 
Great Observatories Origins Deep Survey (GOODS), to date the deepes t wide-area mul 



tiwavelength survey carried out with Spitzer, Hubble, and Chandra (jDickinson et 



2003; iGiavalisco et aLl 12004b . By targeting the Spitzer Legacy, and later the Hubble 



Treasury, observations on the pre-existing Chandra deep fields, we leveraged substantial 
investments of observing time, probed AGN demographics at the peak of quasar activity 
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[z ~ 1 — 2), and enabled a wide range of other science, described in the Astrophysical 
Journal Letters special edition of January 10, 20040 

The GOODS data are deep enough to detect AGN to very high redshift (z £ 6), and 
the volume sampled ensures sizable AGN samples out to z ~ 3. To sample larger volumes 
and to search effectively for rare objects like AGN or massive galaxies at high redshift, 
in 2002 we designed the MUSYC survey (Multiwavelength Survey by Yale and Chile), 
covering one square degree in two equatorial and two southern fields. Three of the four 
regions had already been observed extensively, including the Extended Chandra Deep 
Field South (ECDF-S) and the Extended Hubble Deep Field South (EHDF-S), thus 
leveraging substantial investments with Chandra, XMM-Newton, HST, and the largest 
grou nd-based telescopes|j Soon after starting MUSYC we helped start the COSMOS sur- 
vey ( Scoville et aZ.ll20076h . a 2-square degree field centered at 1000+02 that has now been 
imaged extensively with HST, Spitzer, XMM, and Chandrafij The unprecedented combi- 
nation of area and depth allow a wide variety of science, described in the September 2007 
special edition of the Astrophysical Journal Supplement. O ther releva nt multiwave length 
surveys include the Lock man Ho le (Hasi nger et "Zl l200ll). CLASXS |Yang et aZj l2004lh 
the Extended Groth Strip (jPavis et al.ll2007l). SWIRE dLonsdale et aZ.ll2003h. XBOOTES 



( Hickox et al.ll2006l). HE LLAS2XMM fealdi et al\\200±. SEXSI dHarrison et aLll2003l). 
CYDER (jTreister et a/.l[2005h . CHAMP (|Kim et a/.ll2004f) . and AMSS ()Akivama et al 
20031) . 



3. Finding Obscured AGN 

The question we set out to answer with GOODS, MUSYC, and COSMOS was, "Is 
there a substantial population of obscured AGN that is missed by traditional optical 
surveys?" To answer this question we need to sample the AGN population at z ~ 1 — 2, 
at the peak of the number density. The GOODS survey, an d the Chandra Deep Field 
South (CDF-S) X-ray survey on which GOODS piggy-backed (|Giacconi et were 



designed to sample the AGN population at z ~ 0.5 — 2, which includes the AGN that 
make up the X-ray background. Luminous AGN at higher redshifts could certainly be 
detected but the volume surveyed is too small to expect to see a reasonable number of 
them. 

Early results in the CDFS and oth er deep X-ray fields showed there was a populati on 
of optically faint har d X-ray sources ( Alexander et al. 2001 : Franceschini et al. 2002 6l ) ' 



Mai nieri et al. I l2005h . which collectively comprised a large fraction of the integrated X- 



ray background intensity. However, as optical counterparts were identified and spectra 
obtained, several apparent problems emerged. First, the redshift distribution peaked at 
relatively low values, z ^ 1, lower than expected from the early population synthesis 



t For more details see http://www.stsci.edu/science/goods 
$ The four MUSYC fields include the E(J DF-S, tor which substantial additional Chandra ob- 
servations were later acquired by N. Brandt (lLehmer et aZ.|[2005l : IVirani et ql.ll2006l ); the EHD- 
F-S, for which there not yet any X-ray data; a field centered on a z — 6.3 quasar at 1030+05, 
which has relatively deep XMM data; and a new field centered at 1256+01 ("Castander's Win- 
dow"), an equatorial region with low 100-micron dust emission. The first two have now been 
imaged with Spitzer (IRAC and MIPS) and HST (ACS and NICMOS), and all fi elds are acces- 
sible w it h ALMA. For more details see http://www .a stro.yale.ed u/MUSY C an d lGawiser et all 
( 200661) : lQuadri et all (|2007l h iGawiser et al\ lf2006d ): lvan Dokkum et al\ (|2006D . 

% COSMOS is now one of the largest and best covered deep multiwav elength survey fields , 
representing the investment of thousands of hours of major telescope time (|Scoville et al. 2007a; 
ISanders et aZll2007l ; lHasinger efalll2007i : ICapak et aZll2007l ; iLillv et aZj | 2007f ). For more details, 
see http:/ /cosmos. astro. caltech.edu 
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models for the X-ray background. Second, the fraction of X-ray sources that were identi- 
fied as obscured, either because of high Nh or absence of broad lines, was less than the 
canonical 3/4 seen at low redshift, and appeared to decline with redsh ift rather than in - 
crease, as specified by the best population synthesis mo del at that time dGilli et al. 2001 ) 



A number of authors pointed o ut th ese problems (e.g..lMainieri et aLll2005t IRosati et al 
2002; iBrandt fe Hasingerll2005f ). and lFranceschini et all (|2002qI ) suggested that "the uni- 



fication scheme based on a simple orientation effect fails at high redshifts" and that the 
production of the X-ray background by a collection of obscured AGN "requires major 
revision." 

At th e same time, these faint re d X-ray AGN were copious emitters of infrared ra- 
diation ( Treister et all [20041 12006 ). fully consistent with the unification paradigm. It 
also became apparent early on that, as for optically-selected quasars, the evolution of 
X-ray selected AGN is luminosity dependent, with high-luminosity AGN evolving e arlier 
and more rapidly than low-luminosity AGN ( Ueda et all 120031 : ICowie et al. 2003f ): this 
luminosity-dependent density evolution had not been incorporated in earlier population 
synthesis models for the X-ray background. Finally, we suspected that selection effects 
could play a significant role in affecting the redshifts and optical identifications of survey 
sources. 

Accordingly, we developed a comprehensive quantitative approach, based on the uni- 
fication scenario and incorporating the most recent, best luminosity function and evo- 
lution, to predict the number counts and redshift distributions at any wavelength (for 
the moment, infrared through X-ray, though it could be generalized) for surveys of ar- 
bitrary area, depth, and wavelength. Here we describe the quantitative interpretation of 
the multiwavelength data from GOODS, MUSYC, COSMOS and other deep multiwave- 
length surveys, including the dependence of the obscured fraction of AGN on luminosity 
and redshift. Specifically, we show that most AGN are obscured at all redshifts; that the 
fraction of obscured AGN decreases with luminosity; and that it increases with redshift, 
at least out to redshift z ~ 2.5, an epoch of substantial black hole and galaxy growth. 
Deep infrared and hard X-ray surveys over larger areas will be needed to reach higher 
redshifts and to probe fully the co-evolution of galaxies and black holes. 

3.1. Connecting X-Ray, Optical, and Infrared Surveys 

Our approach was to connect surveys at different wavelengths by assuming something sen- 
sible about AGN spectral energy distributions (SEDs), then combining those with well- 
measured luminosity functions and evolution to understand the sourc e counts and redshift 
distributions of AGN selecte d at a given flux limit at any wavelength ( Treister et a/.ll2004 . 



2006; iTreister fc Urrvl 120061: an alogous approaches have been taken by iBallantvne et 



20061 and iDwellv fe Page 2006). Simu ltaneously, w e constrain the same AGN population 
to fit the X-ray background ([Treister fe Urrv 2005f ). An alternative approach is to model 
only the X-ray spectra of AGN a nd to fit the X-ray background alone with a mixture of 
obscured and unobscured AGN ( Comastri et al. Il995l : iGilli et al. l200ll 120071 ) ; this con- 



strains the demographics but does not connect the X-ray sources to those detected at 
optical and infrared wavelengths — and thus does not use those additional constraints on 
the AGN demographics, nor does it allow a quantitative estimate of the important effect 
of optical or infrared flux limits on the survey content or spectroscopic identifications. 

Briefly, our procedure was as follows: We started with the underlying AGN demograph- 
ics, described by an AGN luminosity function that incorporates dependence on absorb- 



ing c olumn density, and the luminosity-dependent evolution of this function (jUeda et al 



l2003h . Because this luminosity function is based on hard X-ray observations, it is rela- 
tively free of bias against obscured AGN. 
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Figure 2. The infrared spectra of AGN were modeled by dust emission (Nenkova et aZ.|[2002l ') 
from a clumpy torus geometry (left) that, assuming random orientations and adjusting the 
torus geometry to give a 3:1 ratio of column densities above:below Nh = 10 22 cm - , yields an 
Nh distribution consistent with various observational estimates (see text for details). (Right) 
Two examples of model SEDs (lines), which are fully determined from Lx and Nh- These fit 
very well the observed SEDs of local unobscured (datapoints, top) and obscured (bottom) AGN, 
with no free parameters. The composite model SEDs include infrared dust emission, reddened 
quasar spectra (keyed to Nh value) , and an L, host galaxy, linked to the X-rays by the known 
optical-to- X-ray ratio (which depends on Lx)- 



We then developed a set of SEDs, based on the unification paradigm, that represent 
AGN with a wide range in intrinsic luminosity and absorbing column density (param- 
eterized in terms of the neutral hydrogen column density, Nh, along the line of sight). 
Specifically, at optical wavelengt hs (A=0.1-l microns), we used a Sloan Digital Sky Sur- 
vey composite quasar spectrum ([Vanden Berk et a/.ll200l[ ) plus Milky- Way-type redden- 
ing laws and a standard dust-to-gas ratio to convert Nh to Ay; we also added an L» 
elliptical host galaxy, which is the dominant component for heavily obscured AGN. In 
the X-ray (E >0.5 keV), we assumed a power-law spectrum with photon index T = 1.9, 
typical of unobscured AGN, absorbed by column densities in the range log Nh = 20- 
24 cm' 2 . To describe the infrared part of the SEDs (A >1 micron), which in the unifica- 
tion paradigm includes radiation from dust hea ted by absorbed ultra violet through soft 
X-ray photons, we used dust emission models bv lNenkova et al. (2002) in a clumpy torus 
geometry (Fig. [2^), converting to Njj from viewing angle assuming random orientations. 
The re sulting Nh distribution is completely consistent with various observational esti- 
mates (lUeda et all 120031: iDwellv fc Pagejl2006t [Risaliti et a/.|[l99i IComastri et a/.|[l99l 
Tozzi et aLll2006l : iGilli et a/.ll2007i) . AGN models with the same intrinsic X-ray luminosi- 



ties were normalized at 100 microns. To connect the ultraviolet and X-ray parts of the 
SED we used the standa rd dependence of X-ray to optical luminosity ratio on luminosity 
(e.g.. lSteffen et a/.ll2006l ). 

The SED models, parameterized in terms of Lx (as a proxy for intrinsic luminosity) 
and N H (which depends only on viewing angle and torus geometry), describe extremely 
well the local population of AGN. There is some freedom in the choice of torus geometry, 
of course; we selected an aspect ratio that would produce three times as many AGN with 
Nh greater than 10 22 cm -2 co mpared to smaller column densities (Fig. [2^,), which is 
roughly the observed local ratio ( Risaliti et al. 19991 ). Figured shows the observed SEDs 
of two local AGN, one unobscured (top) and one heavily obscured (bottom), compared 
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Hard X-ray [2-8 keV] Band Flux z Band Magnitude 

Figure 3. X-ray (left) and optical (right) counts observed in the GOODS fields (data points) 
agree well with a simple un ification model (solid line), even with no dependence of torus geometry 
on luminosity or redshift (Tr eister et qZj|2004 ). Inclusion of the luminosity dependence (found 
later in larger AGN samples) gives essentially the same result for the GOODS survey. Unobscured 
AGN (Nh < 10 22 cm -2 ) dominate at bright optical magnitudes (dashed line) and obscured AGN 
(Nh > 10 22 cm~ 2 ) at faint z-band magnitudes (dotted line). For R > 24 mag, the approximate 
limit for obtaining decent optical spectral identifications with 8-meter class telescopes, the vast 
majority of sources are obscured AGN. 



to the model SED corresponding to the observed Lx and Nh- With no free parameters, 
the fit is remarkably good overall, meaning that this representation accurately reflects 
the distribution of photons across a very wide wavelength range. 

Combining these model SEDs with the AGN luminosity function and evolution yields 
a distribution of sources as a function of spatial location and wavelength. For a given 
wavelength and flux limit, then, our calculation produces the expected number counts 
and redshift distribution of AGN. This population can be filtered by multiple flux limits 
(for example, one could look at the X-ray counts and redshift distribution of sources above 
arbitrary flux limits in the X-ray, optical, and infrared bands) and/or by other properties 
(e.g., Nh value). Our procedure was to generate the expected content of real surveys, 
and to compare to observations, in order to constrain the underlying AGN demographics, 
specifically, the AGN luminosity function, its evolution, its dependence on Nh, and the 
family of AGN SEDs (which includes the torus properties). 

3.2. Quantitative Results from the Population Synthesis Model: Number Counts, 
Redshift Distributions, and the X-Ray Background 

In a series of papers we showed that this straightforwa rd model matc h es ver y well the 
observed properties of AGN samples. The first paper. [Treister e£~ ( 2004 ). assumed 



for the sake of simplicity that the torus geometry (and hence the Nh distribution) did 
not depend on luminosity or redshift. Even this simplest population synthesis model fits 
beautifully the observed X-ray and optical counts in the deep GOODS survey (Fig. [3]). 
The source population at faint optical magnitudes, particularly the AGN too faint for 
identification with optical spectroscopy, is dominated by obscured AGN. 

Th e unification-ins pired model described above, like other population synthesis models 
(e.g.. lGilli et a/.ll200ll) . predicts a large population of obscured AGN out to high redshift. 
The predicted peak of the redshift distribution (dashed histogram, Fig. [4]a.) is near z ~ 
1, not much higher than that observed in the GOODS data (shaded/ open histogram, 
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Figure 4. (Left) Observed spectroscopic redshift distribution (hatched histogram) and pho- 
tometric redshift distribution (open histogram) for the GOODS-North field, compared to that 
expected from the unification model (dashed line). There is poor agreement at high redshifts. 
(Right) Imposing an optical flux limit on the model leads to an expected distribution (dashed 
line) that agrees very well with the observed distribution. 



Fig. but there discrepancy at z > 1. Many of these high-redshift obscured AGN, 
because they are optically faint, are not identified in follow-up optical spectra of X-ray 
sources. This can be clearly seen imposing the spectroscopic limit, R ~ 24 mag, on 
the expected redshift distribution (dashed histogram, Fig. [4})); in this case, the observed 
and expected distributions agree well because the higher-redshift AGN remain largely 
unidentified. Equivalently, AGN without spectroscopic identifications will preferentially 
be faint, obscured sources. 

The advent of multiple large AGN samples with high spectroscopic completeness al- 
lowed B arger et al. (2005) to deduce a dependence of obscured fraction on optical lumi- 
nosity. We incorporated a simple functional form of that dependence in our model — 
namely, a linear transition between 100% of AGN with Lx = 10 42 ergs/s being obscured 
and no AGN with Lx — 10 46 ergs/s being obscured. Taking into account the selection 
effects due to flux limits, the model then matches the observed dependence very well 
(Fig. [TTj) . The resulting number counts and redshift distributions for GOODS do not 
change, since only a restricted lumino sity range (primarily 10 43-44 ergs/s) is probed by 
that survey. In all work subsequent to iTreister et all (|2004f ) — on infrared counts, X-ray 
background, and evolution of obscuration (see below) — we incorporated this luminosity 
dependence. 

Our population synthesis model predicts, with very little freedom, the spectrum of the 
X-ray background. We simply sum the X-ray emission from all the AGN. There is some 
freedom in that we (like everyone else who does this kind of calculation) have to model 
the X-ray spectrum. We used an absorbed power law plus an iron line and a Compton- 
reflection hump, consistent with the AGN spectra that are well measured locally; we also 
assumed solar abundances and used the Njj distribution described above for log Njj = 20- 
24 cm" 2 . The Nh distribution for higher column densities (Njj > 10 24 cm" 2 ) is poorly 
constrained at present, so we extrapolated from Njj < 10 24 cm -2 with a flat slope (again, 
much like what others have assumed). Beyond this, one can adjust the metallicity (this 
increases the hard X-rays, holding everythi ng else constant) or the include a small range 
of power- law slopes (ditto: iGilli et aZj|2007t ). but we did not feel the data could constrain 
those parameters and so left them fixed. In any case, the integral constraint of the X-ray 
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Figure 5. The predicted X-ray background spectrum (heavy line; iTreister fc Urry||200"5T ) from 
the unification- inspired population synthesis model agrees very well with the observed data (data 
points). (Left) Obscured AGN (thin, red, peaked curve) dominate at high energies (E > 10 keV), 
while unobscured AGN (thin, blue, broad curve) dominate in the Chandra and XMM bands 
(2 keV < E < 10 keV). There is little freedom in the fit parameters, especially at low energies, 
so this agreement is remarkably good. (Emission at E > 100 keV is dominated by blazars, which 
are not modeled here.) The X-ray background is dominated by moderate luminosity AGN, in the 
range iq 43-44 ergs/s (middle panel) and by intermediate redshifts, z ~ 0.5 — 1.5 (right panel). 



background is an excellent check on whether the demographics of AGN assumed here 
is realistic. Figure [5] shows that, with almost no free parameters, the data are very well 
fit indeed. (The normalization of the background at E > 10 keV is discussed below in 
§ 3.3.) 

Another strong prediction of the Treister model concerns the infrared data, which at the 
time the model was developed had not yet been taken. We calculated the expected Spitzer 
counts in IRA C and MIPS 24-micr on bands; Figure [5] shows the excellent agreement with 
observations ( Treister et al. 20061 ). Small discrepancies at the faint end are due to the 



overly simple assumption of a single host galaxy magnitude (which dominates at low 
fluxes), rather than a distribution. 

From the Spitzer observations, we calculate the minimum AGN contribution to the ex- 
tragalactic infrared background, obtaining a lower value than pre viously estimated, rang- 
ing from 2% to 10% in the 3-24 micron range ( Treister et aLl[2006h . Accounting for heavily 
obscured AGN that, according to our population synthesis model, are not detected in 
X-rays, the AGN contribution to the infrared background increases b y ~45%, to ~3-15% . 
FigureEJi shows the AGN contributions deduced from GOODS data (|Treister et aZ.ll2006l) 
compared to other estimates and to the (uncertain) total extragalactic background light, 
as a function of infrared wavelength. The GOODS measurements place the strongest 
constraints to date on the AGN contribution to the extragalactic background light, in- 
dicating that stars dominate completely over AGN at infrared wavelengths. 

The fraction of sources that are AGN rises sharply with 24 /im infrared flux (Fig.[7)3). 
Thus in deep surveys like GOODS or COSMOS, AGN are a small fract ion of the total 
infrar ed source population, while in high flux limit surveys like SWIRE ( Lonsdale et all 
20031) . they constitute a much higher fraction. AGN detected in large-area, shallow sur- 



veys are on average closer and/or more luminous than AGN found in deep pencil-beam 
surveys. We show in Figure [TJd that this very strong dependence of AGN fraction on 
infrared flux is not an artifact of the X-ray flux limit. Specifically, using our AGN popu- 
lation synthesis model, we plot the number of AGN as a function of infrared flux including 
those too faint to be detected in X-rays in the Chandra Deep Fields (open circles). Clearly 
the trend is independent of X-ray flux limit. 
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Figure 6. (Filled circles) Observed cumulative infrared flux distributions for X-ray detected 
AGN in the GOODS-North and South fields in three IRAC bands, and for the North field 
only in the MIPS 24-micron band. Error b ars correspond to the 84% confidence level in the 
number of sources in that bin (Gchrels 1986). Because this is a cumulative plot, the bins are not 
independent. (Solid line) Power- law fit to the observed number counts. (Dashed line) Expected 
infrared flux distribution (Trei ster et al.\\200(i ) based on our population synthesis model, taking 
into account the X-ray flux limit in the GOODS fields and the luminosity dependence of the 
obscured AGN fraction, and correcting for the effects of undersampling at the bright end caused 
by the small volume of the GOODS fields. (Dotted line) Expected infrared flux distribution not 
considering the X-ray flux limit (i.e., all sources expected from extrapolating the AGN luminosity 
function). In general, the model explains the normalization and shape of the counts quite well, 
especially considering the poor statistics at the bright end. At the faint end, where the host 
galaxy light dominates, the model diverges more significantly, since the model assumptions (no 
distribution in host galaxy magnitude or spectrum) are clearly far too simple. 



3.3. Integral and Swift Hard X-Ray Surveys 

If the population synthesis model presented here is correct, ^50% of AGN are currently 
missed even in deep X-ray surveys with Chandra or XMM. These are very obscured 
AGN, many of them Compton thick (i.e., Njj > 10 24 cm -2 ), so all but the hardest X-ray 
surveys are biased against them. Fortunately, hard X-ray instruments on the INTE- 
GRAL and Swift satellites can now reach fluxes below ~10 -11 ergs cm~ 2 s" 1 at energies 
above 20 keV. Serendipitous AGN surve ys at these energies cove ring almost the full sky 
have been carried out using Swift /BAT dMarkwardt et aLl l2005h and INTEGRAL/IBIS 
( Beckmann et al. 2006HSazonov et a/1l2007 ). yielding samples of ~100 AGN each. These 
surveys provide an unbiased view of the AGN population, independent of the amount of 
obscuration, although at present the sensitivity is sufficient to reach only local popula- 
tions. 

Figure^ shows the logN-logS distribution from the INTEGRAL catalog (Sa zonov et al 



120071) of AGN confirmed as Compton thick with Nh measurements fro m X-ray data 
Clearly the original population synthesis model of Treister fc Urrvl ( 20051 ) overpredicted 
the density of Compton-thick AGN by a factor of ~4. This is not too surprising, as 
the model included assumptions that, at the time of the publication of that work, were 
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Figure 7. (Left) Extragalactic infrared background intensity as a fu nction of wave- 
length . (Shaded region) Allowed values for t he total intensity compi led by lHauser fc D wck 
(|200ll ). 1 -4 y,m: Mea s ureme nts reported by iMatsumoto et al\ (12005 ) . 4- 8 /im: Lower lim- 
its from iFazio et al.\ (|2004h . upper limits from iKashlinskv et all ( 19961 1 a t 4.5 tim and 
IStanev fc Franceschinil (j 1998f t at 5.8 and 8 ^m. Measurement at 24 /im from iPapovich et al.\ 
(|2004T ). (Solid circles) Integrated infrared emission from X-ray-detected AGN in the GOODS 
fields, with la error b ars from sourc e number statistics (which dominate over measured flux 
errors) calculated as in iGehrelsl (|T986). (Triangles) Integrated AGN emission from the models of 
Trei ster fc Urrvl (|2005l ) if only X-ray-detected AGN are considered (lower) or including all AGN 
(upper). (Dotted lines) Expected AGN integrated emission from the AGN models of Silva et al. 
(2004; upper line) and Xu et al. (2001; lower line). (Open squares) Expect ed AGN integrated 
emission at 15 /im from the extrapolation to fainter fluxes of IRAS and ISO (jMatute et qZ.|[2002l 
[M021 ; corrected by a f actor of 5 to include the contribution from obscured AGN) and ISO 
only ()Fadda et qZ.ll2002l [F02]) observations. (Right) The fraction of sources that are AGN rises 
sharply with 24 (im infrare d flux (filled c ircles). Vertical error bars show the la Poissonian errors 
on the number of sources (Gchrcls 1986). Also shown is the contribution corrected by the AGN 
expected to be missed by X-ray selection (open circles), as estimated using our AGN population 
synthesis model; this shows the same strong dependence on infrared flux, indicating that that 
dependence is not a selection effect induced by the X-ray flux limit. 



unconstrained — for example, the reflection fractioijj] was assumed to be unity, and the 
number of Compton-thick AGN was simply extrapolated with a flat slope from the Nh 
distribution at lower column densities. 

Now we explore the constraints on these quantities that come from the new Swift 
and INTEGRAL surveys. First, we define a "Compton-thick AGN correction factor," 
which simply multiplies the original flat-extrapolation assumption to match the observed 
number density of AGN with Njj > 10 24 cm -2 . According to a x 2 minimization, the 
best-fit value for the Compton-thick AGN factor is 0.25; this produces a good agreement 
between model and observations, with a reduced \ 2 °f 0-3- 

Second, we consider the ratio of direct and reflected X-rays. The spectrum of Compton- 
thick AGN at h igh energies is dominated by the Compt on reflection component (e.g ., 
Matt et aZ.ll2000h . which has a strong peak at E ~30 keV ([Magdziarz fc Zdziarskilll995l ). 



The observed spectrum of the X-ray background, which we now understand as the inte- 
grated emission from previously unresolved AGN, also has peak at about the same energy 
(|Gruber et aHll999h . The normalization of the reflection component relative to the direct 



f The reflection fraction is the geometrical factor describing the solid angle, relative to 2n, 
subtended by cold reflecting material as seen from the primary X-ray source. 
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Figure 8. (Left) Lo gN-logS relation for Co mpton-thick AGN only. The dashed line shows the 
original prediction of lTreister fc Urrvl (|2005f) . which assumed a reflection component normaliza- 
tion of 1 and a Compton-thick AGN factor of 1 (i.e., as many Compton-thick AGN as Comp- 
ton-thin AGN in the next lowest decade of the Nh distribution). Data points are from the 
INTEGRAL survey of [S azonov et al.\ (|2007l ). The population synthesis model can be brought to 
agreement with the observations if there are a factor of 4 fewer Compton-thick AGN (i.e., the 
Compton-thick AGN factor is 0.25). (Right) Compton-thick AGN factor versus normalization 
of the Compton reflection component. The gray region shows the values of these parameters 
tha t produce a space de nsity of Compton-thick AGN consistent with the observed value in 
the ISazonov et al\ (|2007f ) sample, considering 1-a statistical fluctuations. The allowed values 
of the parameters needed to fit the intensity of the X-ray background at 20-40 keV, assum- 
ing a 5% uncertainty in each case, are shown by the blue region for the original iGruber et al.\ 
(|1999|) measurements, the red region for the INTEGRAL values and by the green region for 
the | Gruber et al\ (|1999D points increased by 40% (which was the assumption in ITreister fc Urrvl 
2005). A large value of the reflection component is required to obtain values consistent with the 
latter X-ray background intensity, inconsistent with observations of the re flection compo n ent in 
individual AGN. Both the reflection value and the renormalization of the IGruber et al\ l| 19991 ) 
X-ray background intensity must be lower. 



emission is known only for a few nearby AGN, mostly from BeppoSAX observations (e.g., 
Perola et al. 20021 ). and therefore this parameter is usually taken as an assumption in 



AGN population synthesis models that can explain the spectral shape and intensity of 
the X-ray background. The resulting peak X-ray background intensity depends on both 
the assumed space density of Compton-thick AGN and the normalization of the reflection 
component, R; while satisfying the overall intensity constraint, one can trade increased 
reflection for fewer Compton-thick AGN, or vice versa. Figure [5}d shows the allowed re- 
gions in terms of the Compton-thick AGN factor versus reflection parameter. That is, 
the shaded regions show the values of these two parameters that produce an integrated 
X-ray background intensity in the 20-40 keV range, the region most affecte d by both the 



numb er of Compton-thick AGN and R, consistent with: (green region) the IGruber et 



(1999) data increased by 40% (the maximum suggest ed renormalization to match higher 
estimates with imaging instruments at lower energies; iTreister fc Urry||2005i ): (red region) 
the recent INTEGRAL X-ray background intens ity measured using E arth occultation by 
Churazov et al. ( 20061 ): (blue region) the original Gruber et al. ( 1999f ) measurements, not 



renormalized. In each case uncertainties in the X-ray background intensity were assumed 
to be 5%. The gray shaded region shows the parameter space that produces a density of 
Compton-thic k AGN consistent with the observations in Fig. [5^. 

Previously. ITreister fc Urry ( 2005 ) assumed an X-ray background intensity consistent 
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Figure 9. (Left) Area versus hard X-ray flux for each of the seven surveys combined in this work 
(light colored lines) and for the total sample (thick black line). (Right) Total effective area as a 
function of X-ray flux and optical magnitude (for R—24, 26 and 28 mag), taking into account 
the spectroscopic incompleteness of each survey (see text for details). These curves are used 
to compute the expected fraction of identified obscured AGN for an intrinsically non-evolving 
ratio. 



with the Gruber et all ( 1999f ) value increased by 40%, which was well fit by a Compton- 
thick AGN factor of 1 (flat extrapolation of the Nh distribution; Fig. [5]). However, 
such a high value of the Compton-thick AGN factor is clearly inconsistent with the new 
observational constraints on the density of Compton-thick AGN. So, either the intensity 
of the X-ray backgro und is lower, as suggested by recent INTEGRAL measurements 
( Churazov et a/.ll2006l ). or the average value of the reflection parameter is high, R ~2, or 
some combination of the two. Observations of individual sources seem to indicate that 
such a high value for the reflection compone nt is unlikely. From a sample of 22 Seyfert 



galaxies, excluding Compton-thick sources, iMalizia et al\ (|2003l ) concluded that both 
obscured and unobscured sources have similar reflections component with normal i zation 
values in the 0.6-1 range. A similar value of R ~1 was reported by lPerola et al. (2002) 
based on BeppoSAX observations of a sample of nine Seyfert 1 galaxies. Although with 
large scatter, normalizations for the average reflection component of 0.9 for Seyfert 1 and 
1.5 for Seyfert 2 were measu red from BeppoSAX observations of a sample of 36 sources 
(jDeluit fc Courvoisierl 120031 ). Therefore, a value of R ~1 for the normalization of the 
reflection component, required by both the observed Compton-thick AGN space density 
and the X-ray background intensity reported by INTEGRAL and HEAO-1, is in good 
agreement with the observed values in nearby Seyfert galaxies. 



4. Evolution of the Obscured Fraction of AGN 

The X-ray background, being an integral constraint, is not a strong probe of the 
fraction of AGN that is obscured (as we showed in the previous section), much less of the 
evolution of that fraction with cosmic epoch. As shown in Figure O a simple population 
synthesis model in which obscured and unobscured AGN have the same evolution (or 
equivalently, the fraction of AGN that is obscured does not evolve) is fully consistent 
with the spectrum of the X-ray background. Now, however, the large X-ray samples 
that have become available in the past few years, spanning a range of depths and with 
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Figure 10. Fraction of optically identified AGN as a function of optical magnitude (data points), 
for four of the seven surveys comprising our super-sample. All seven are well described by a 
simple linear increase to a constant fraction at bright magnitudes (lines). This allows us to 
derive the effective survey area as a function of both X-ray flux and optical magnitude. 

high spectroscopic completeness, all ow us to determine wh ether and how the fraction of 
obscured AGN depends on redshift ( Treister fc Urry 20061 ). 

To study the evolution of the obscured AGN fraction, one needs to distinguish between 
the effects of redshift and luminosity, which are c orrelated in any flu x-limited sample, 
Wide area, shallow X-ray surveys (e.g., XBOOTES: fHickox et aLll2006h sample moderate 
luminosity AGN at low redshifts and only hi gh luminosity sources up to high redshifts, 
while deep pencil-beam surveys (e.g., CDFS; Giacconi et a/.|[2002T) find moderate lumi- 
nosity AGN out to high redshifts but lack rare, high-luminosity sources because of the 
small volume sampled. Combining the two extremes covers the luminosity-redshift plane 
effectively. 

We generate an AGN su per-sample comprisin g seven l arge surveys with h igh iden- 



tification fractions: AMSS (Akiyama et al. 



LAS2XMM (jBaldi et all l2002h . CLASXS ilVanf iTu~i\ 1200 J . CYDER. fTreister et 



2003), SEXSI (H 



arrison et all 2003), HEL- 
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2005) , and the two Chandra deep fields (jGiacconi et al. 2002t lAlexander et al. 2003 ). 
This super-sample contains a total of 2341 hard X-ray-selected AGN (X- ray sources with 
Lx > 10 42 ergs/s), the largest such sample to date by a factor of ~4 (jTreister fc Urry 

2006) . It spans a range of luminosities at each redshift, over a broad redshift range. The 



total area of this super-sample as a function of X-ray flux is shown in Figure [9^. 

More than half the super-sample is optically identified. We define an AGN as unob- 
scured when there is evidence for broad emission lines in their spectra, and as obscured 
AGN otherwise. The "obscured fraction" is then the number of obscured AGN divided 
by the total number of AGN. For obvious reasons, the identified fraction of any survey 
depends on the brightness of the optical counterparts. We parameterized the identified 
fraction of each of the seven surveys with a simple function t hat is constant at bri ght 
magnitudes and falls linearly to faint magnitudes (Fig.[l0l see lTreister fc Urry 20061 for 
details). We then weight the area versus X-ray flux curve (Fig. [9£i) by the complete- 
ness at each optical magnitude, for each survey, and sum those to get the effective area 
of the super-sample as a function of both X-ray flux and optical magnitude (Fig. [SJd). 
This allows us to calculate the expected numbers of optically identified AGN for the 
super-sample, i.e., we can now correct for both X-ray and optical spectroscopic limits. 

As discussed ear lier (§ 13. 2[) . the fraction of obscured AGN depends on luminosity 
( Barger et aLll2005l) . Our population synthesis model assumed a linear transition between 
100% obscured fraction at Lx = 10 42 ergs/s and obscured fraction at Lx = 10 46 ergs/s. 
Taking into account the selection effects due to flux limits (Fig. [9]d) , this assumption 
(black line, Fig. Illj) matches the observed dependence very well in the present super- 
sam ple (black line, Fig. Illj) . The somewhat shallower luminosity dependence adopted 
by ( Gilli et aZ.ll2007t blue dashed line is their assumed dependence, red dotted line incor- 
porates flux limits), in contrast, does not describe the obscured fraction well for high- 
luminosity AGN, probably because their model was constrained by the X-ray background, 
which is dominated by moderate luminosity AGN. 

Our population synthesis model fits the X-ray, optical, and infrared counts of AGN; 
the X-ray background (modulo the trade-off between the number of Compton-thick AGN 
and the reflection fraction of each); the hard X-ray counts measured with Swift and 
INTEGRAL (ditto); and the observed luminosity dependence of the obscured fraction of 
AGN in the largest AGN sample to date. What can it tell us about the evolution of this 
obscured fraction? 

Figure [T2"a shows the observed fraction as a function of redshift (black data points). As 
many have noted previously, the observed fraction declines with redshift, from roughly 
3/4 locally to ~ 1/3 at z ~ 4. This has been interpreted to mean that obscured AGN 
are rare at high redshift. However, the lines in Figure [T2a show the expected decline 
for an underlying population whose obscured fraction is actually constant with redshift, 
calculated for our super-sample using the appropriate corrections for X-ray and optical 
flux limits and spectroscopic completeness. That is, an even steeper decline is expected 
in the observed samples even when the underlying population does not change at all with 
redshift. 

Three significant selection effects cause this strong decline. First, obscured AGN have 
smaller X-ray fluxes, so there is a bias against their detection in the X-ray sample in 
the first place. This is a small effect, particularly because it becomes relatively less 
important with increasing redshift (for which rest-frame emission is in an increasingly 
harder X-ray band affected only by higher and higher column densities). Second, obscured 
AGN are fainter in the optica l, so there is a bias against spectroscopic identifications 
for z ^ 1 ( Treister et al. 20041 ): this is not important for samples with highly complete 



spectroscopic identifications. Third and quite important, the obscured fraction depends 
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Figure 11. Fraction of obscured AGN in the super-sample as a function of hard X-ray lumi- 
nosity. (Black data points) Observed distribution for the super-sample described in this work. 
(Red data points ) Compilation of G. Hasinger (2007, in prep.) used in the work reported by 
Gilli et al\ (|2007l ). (Black solid line) Expected luminosity dependence for the Treister popula- 
tion synthesis models, taking into a ccount the X-ray flux l imit and spectroscopic completeness 
of the super-sample described here (jTreister fc U rry 2006). (Blue dashed and red dotted lines) 
Intrinsic and bias-co r rected (i.e., expected) luminosity dependences for the population synthesis 
model of Gill i et all (|2007T ) . This assumed lum inosity dependen ce is not a good description of 
AGN demographics at higher luminosities; the iGilli et aZ.I [20071 model is constrained primarily 
by the X-ray background intensity, which is dominated by moderate luminosity AGN. 



strongly on luminosity (Fig. [TTjl . Given the luminosity- redshift correlation inherent in 
flux-limited samples, the mean AGN luminosity increases with redshift and therefore the 
obscured fraction that is observed decreases — even if the same population of obscured 
lower-luminosity AGN is present. Our analysis shows that this is the dominant selection 
effect. It is an important effect even in AGN samples with 100% complete spectroscopic 
identifications, as indicated by the blue line in Figure[T2k. Simply put, high redshift AGN 
samples are biased to higher luminosity, and thus contain lower obscured fractions, even 
if there is no underlying evolution of the ratio between obscured and unobscured AGN 
at all. 

The observed more gradual decline in obscured fraction of AGN actually implies an 
increase in obscured fraction with redshift. Figure [12b shows the data relative to (i.e., 
divided by) the expectation for a non-evolving population, for the super-sample and the 
two sub-samples with higher identification fractions. The effect is largely independent 
of the completeness of the optical identifications. Fitting the increasing fraction with a 
simple power-law dependence on redshift, oc (1 + z) a , gives a good fit for a = 0.4 ± 0.1. 
This means that the fraction of AGN at redshift z ~ 4 that are obscured is observed to be 
twice as high as would be the case were the intrinsic fraction constant. AGN obscuration 
is substantially greater in the young Universe. 



1G 
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Figure 12. (Left) Observed fraction of obscured AGN as a function of redshift. (Data points) 
Measured fraction for our super-sample of 1229 optically-classified AGN (black points), and 
for sub-samples cut at brighter optical magnitudes to increase the spectroscopic completeness 
to 73% (red points) and 93% (blue points). (Solid lines) Expected fraction as a function of 
redshift for an intrinsically non-evolving ratio, taking into account the effects of spectroscopic 
incompleteness and flux limits for the super-sample and two sub-samples. (Right) Fraction of 
obscured AGN relative to the expectations for a non-evolving obscured AGN ratio (dotted line), 
incorporating the effects of spectroscopic incompleteness. A significant increase with redshift is 
clearly seen, independent of the spectroscopic completeness of the sample. Dashed lines show an 
intrinsic evolution of the form (l+z) a , for a = 0.4 ± 0.1. 



5. Summary 

GOODS, MUSYC, and other deep multiwavelength surveys provide overwhelming ev- 
idence for a large population of obscured AGN that dominate AGN demographics out 
to high rcdshifts. Optical surveys are biased against detecting these objects, and even 
hard X-ray surveys, which are considerably less biased, suffer strong selection effects, 
primarily due to the luminosity-dependence of obscuration. 

Taking these effects into account quantitatively, with a realistic, well-constrained pop- 
ulation synthesis model that uses AGN spectral energy distributions based on a unifica- 
tion paradigm, we deduce that the ratio of obscured (Nh > 10 22 cm~ 2 ) to unobscured 
(Nh < 10 22 cm~ 2 ) AGN is roughly 3:1 locally (integrated over all luminosities), and 
increases with redshift. Low-luminosity AGN (10 42 ergs/s < Lx < 10 44 ergs/s) are much 
more likely to be obscured than high-luminosity AGN (Lx > 10 44 ergs/s). 

To the extent that our assumed infrared through X-ray spectral energy distributions 
are reasonable, and our assumed Njj d istribution is rea sonable (it is essentially the same 
as that used by others in the field; see iGilli et al. 2007 for a comparison of the different 
distributions), these results are completely robust. 

How might one avoid the biases inherent in optical and X-ray surveys? In principle, 
far-infrared surveys are unbiased because the absorbed energy is re-radiated thermally; 
however, infrared surveys are very inefficient for AGN selection since the infrared sky is 
strongly dominated by starlight (Fig. [7^,). In addition, AGN signatures (such as broad 
emission lines, strong power-law continua, and rapid variability) may well be hidden in 
obscured objects. Thus, identifying complete samples of AGN from far-infrared surveys 
will never be a simple matter, although it can potentially put useful limits on the fraction 
of galaxies with buried AGN. 

Selection effects are very important to take into account for any survey, even those 
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that are 100% identified. For example, consider a deep hard X-ray survey for which all 
sources have optical and/or infrared counterparts with known redshifts and classification. 
Many would assume the survey itself is not strongly biased and that the identifications 
yield the underlying demographics, since no X-ray sources remain unidentified. However, 
missing from the X-ray sample are the most heavily obscured AGN; even more impor- 
tant, unobscurcd AGN are over- represented (relative to their fraction of the underlying 
population), especially at high rcdshift, because of the dependence of obscuration on 
luminosity. This in fact is the dominant effect for existing surveys. 

Therefore, to understand the distribution of black holes in the universe and to estimate 
cosmic accretion rates, the selection effects must be modeled using reasonable assump- 
tions about the underlying population. That in turns yields a picture of the universe 
that matches well our picture of nascent accreting black holes at the centers of dusty, 
star-forming, young galaxies in the early Universe. 

This work would not have been possible without the exquisite data made possible by 
NASA's Great Observatories and by great observatories on the ground. We gratefully 
acknowledge support from NSF grant AST-0407295 and NASA grants NAG5-10301, 
HST-AR-10689.01-A, HST-GO-09822.09-A, HST-GO-09425.13-A, and NNG05GM79G. 
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